#### Results, Figures and Tables Replication
#### The Effect of Reciprocity on Trust: International Cooperation and COVID Aid
#### Andrew Roskos-Ewoldsen, Morgan Ellithorpe, and Brandon Kinne

## Clear the environment
rm(list=ls(all=TRUE))

## Set your working directory
wd <- c("YourDirectoryHere")
setwd(wd)

## Load packages

library(tidyverse)
library(stargazer)
library(ggplot2)
library(viridis)

## Load data

d1 <- read.csv(file.choose())
# You should select the file titled "qualtricssample_final.csv" from the replication folder

### Results, Figures, and Tables from Main Text

## Reproducing Figure 1: Trust by Treatment Group
## Mean Response by Treatment Group Plot 

# Narrow sample to exclude any missing data.
d2 <- d1[complete.cases(d1[ , c('trust_total', 'similaritycond', 'recipcond', 'costcond')]), ] 

# Create labels for treatment groups
d2 <- d2 %>%
  mutate(Reciprocity = factor(recipcond, levels = c(0, 1), labels = c("Reciprocity Violation", "Reciprocity")),
         Costliness = factor(costcond, levels = c(0, 1), labels = c("Low Cost", "High Cost")),
         Similarity = factor(similaritycond, levels = c(0, 1), labels = c("Low Similarity", "High Similarity")))

#Calculate mean responses for each treatment group, as well as confidence intervals
mean_ci <- d2 %>%
  group_by(Reciprocity, Similarity, Costliness) %>%
  summarize(mean_trust_total = mean(trust_total),
            ci_low = mean(trust_total) - qt(0.975, n() - 1) * sd(trust_total) / sqrt(n()),
            ci_high = mean(trust_total) + qt(0.975, n() - 1) * sd(trust_total) / sqrt(n()))

# Create a color palette for plot using Viridis package
pal <- viridis(2, begin=0.2, end=0.8, direction=-1)

# Create the plot using ggplot. 
meanplot <- ggplot(mean_ci, aes(x = Reciprocity, y = mean_trust_total, shape = Reciprocity,
                                colour = Reciprocity)) +
  geom_point(size = 2.5) +
  geom_linerange(aes(ymin=ci_low, ymax=ci_high), linewidth=1, alpha=0.75) +
  facet_grid(. ~ Similarity + Costliness, scales = "free_y") +
  scale_y_continuous(limits = c(0, 10)) +
  scale_color_manual(values=pal) +
  labs(y = "Mean Response") +
  theme_bw() +
  theme(axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        legend.title = element_blank(),
        legend.position = "bottom",
        legend.text = element_text(size = 11)
  )
#ggsave("figure1.pdf", plot = meanplot) #You can save the plot as a pdf using these lines of code.
#dev.off()

## Reproducing Figure 2: The Reciprocity Effect on Trust
## Forest plot for base linear regression model

# Specify what DVs to include
dvs <- c("trust_total")

# Apply lm() function to each DV; use same covariates in each model
ests <- lapply(dvs, function(x) {
  lm(substitute(i ~ similaritycond + costcond + country + age + male +
                  education + vax + familyhadcovid + demdummy + repdummy +
                  polorient + CIfa + MIfa + ISOfa + recipcond, list(i = as.name(x))), data=d1)
})

# Name list elements for DVs
names(ests) <- c("Trust")

# Put in a ggplot-friendly data frame
ests.df <- data.frame(
  betas = unlist(lapply(ests, function(x) x$coefficients)),
  ses = unlist(lapply(ests, function(x) sqrt(diag(vcov(x))))),
  pval = unlist(lapply(ests, function(x) summary(x)$coefficients[,4])),
  var = unlist(lapply(ests, function(x) names(x$coefficients))),
  mod = rep(names(ests), each=length(ests[[1]]$coefficients))
) %>% # Calculate CIs
  mutate(ci_hi = betas + (1.96 * ses),
         ci_lo = betas - (1.96 * ses)
  ) %>% # Flag significant estimates
  mutate(sig = if_else(pval < .05, 1, 0)
  ) %>% #Rename variables to look better in plot
  mutate(var = recode(var,
                      "(Intercept)" = "Intercept",
                      "similaritycond" = "Similarity",
                      "costcond" = "Cost",
                      "country" = "Country",
                      "age" = "Age",
                      "male" = "Male",
                      "education" = "Education",
                      "vax" = "Vaccinated",
                      "familyhadcovid" = "Covid Experience",
                      "demdummy" = "Democrat",
                      "repdummy" = "Republican",
                      "polorient" = "Political Orientation",
                      "CIfa" = "Cooperative Internationalism",
                      "MIfa" = "Militant Internationalism",
                      "ISOfa" = "Isolationism",
                      "recipcond" = "Reciprocity")
  ) %>% #Rearrange variables for more sensible order in plots
  arrange(betas) %>%
  mutate(var = factor(var, levels= rev(c("Intercept", "Reciprocity",
                                         "Similarity", "Cost",
                                         "Cooperative Internationalism", "Militant Internationalism",
                                         "Isolationism", "Country",
                                         "Age", "Male", "Education",
                                         "Democrat", "Republican",
                                         "Political Orientation", "Vaccinated",
                                         "Covid Experience"))))

# Create color palette for plot using viridis package
pal <- viridis(2, begin=0.2, end=0.8, direction=-1)

# Create forest plot
basemodel_gg_v1 <- ests.df %>%
  filter(var != "Intercept") %>% #For excluding certain variables
  ggplot(aes(x=betas, y=var, color=as.factor(sig))) +
  geom_vline(xintercept=0, lty=2) +
  geom_point(size=2, shape=16) +
  geom_linerange(aes(xmin=ci_lo, xmax=ci_hi), linewidth=1, alpha=0.75) +
  scale_color_manual(values=pal) +
  facet_wrap(~mod) +
  xlab("Estimated Coefficients + 95% Confidence Intervals") +
  theme_bw() +
  theme(axis.title.y = element_blank(),
        legend.position="none",
        axis.text.y = element_text(color="black"),
        axis.text.x = element_text(size=8)
  )
#ggsave("figure2.pdf", plot = basemodel_gg_v1) # You can use these two lines of code to save the plot as a pdf
#dev.off()

## Reproducing Figure 3: Reciprocity Conditional on Similarity and Costliness
## Forest plot for interaction models

# Estimate interaction models
ests1 <- lm(trust_total ~ recipcond*similaritycond + country + age +
              male + education + vax + familyhadcovid + demdummy + repdummy +
              polorient + CIfa + MIfa + ISOfa + recipcond, data=d1
)

ests2 <- lm(trust_total ~ recipcond*costcond + country + age +
              male + education + vax + familyhadcovid + demdummy + repdummy +
              polorient + CIfa + MIfa + ISOfa + recipcond, data=d1
)

# Put interaction models in list, rename list elements
ests.inter <- list(ests1, ests2); rm(ests1, ests2)
names(ests.inter) <- c("Similarity", "Costliness")


# Put in ggplot-friendly data frame
ests.inter.df <- data.frame(
  betas = unlist(lapply(ests.inter, function(x) x$coefficients)),
  ses = unlist(lapply(ests.inter, function(x) sqrt(diag(vcov(x))))),
  pval = unlist(lapply(ests.inter, function(x) summary(x)$coefficients[,4])),
  var = unlist(lapply(ests.inter, function(x) names(x$coefficients))),
  mod = rep(names(ests.inter), each=length(ests.inter[[1]]$coefficients))
) %>% # Calculate CIs
  mutate(ci_hi = betas + (1.96 * ses),
         ci_lo = betas - (1.96 * ses)
  ) %>% # Flag significant estimates
  mutate(sig = if_else(pval < .05, 1, 0)
  ) %>% # Rearrange order of mod variables for plot order
  mutate(mod1 = factor(mod, levels=c("Similarity","Costliness")))

# Keep only the variables of interest
kps <- c("recipcond", "similaritycond","costcond", "recipcond:similaritycond",
         "recipcond:costcond")

## Specify color palette
pal <- viridis(2, begin=0.1, end=0.8, direction = -1)

#Create forest plot
gg_inter <- ests.inter.df %>% 
  filter(var %in% kps) %>%
  mutate(var = recode(var,
                      "similaritycond" = "Similarity",
                      "costcond" = "Cost",
                      "recipcond" = "Reciprocity",
                      "recipcond:similaritycond" = "Reciprocity x Similarity",
                      "recipcond:costcond" = "Reciprocity x Cost",)
  ) %>% #Rearrange variables for more sensible order in plots
  arrange(betas) %>%
  mutate(var = factor(var, levels= c("Reciprocity", "Similarity", "Reciprocity x Similarity",
                                     "Cost", "Reciprocity x Cost"))) %>%
  ggplot(aes(x=betas, y=var, color=as.factor(sig), shape=as.factor(sig))) +
  geom_vline(xintercept=0, lty=2) +
  geom_point(size=2, shape=16, position=position_dodge(0.5)) +
  geom_linerange(aes(xmin=ci_lo, xmax=ci_hi), linewidth=1, alpha=0.75, position=position_dodge(0.5)) +
  scale_color_manual(values=pal) +
  xlab("Estimated Coefficients + 95% Confidence Intervals") +
  theme_bw() +
  facet_wrap(~mod1, scales="free_y") +
  theme(axis.title.y = element_blank(),
        legend.position="none",
        legend.title = element_blank(),
        axis.text.y = element_text(color="black"),
        axis.text.x = element_text(size=8)
  )
#ggsave("figure2.pdf", plot = gg_inter, width = 16, height = 10, units = "cm") #you can use these two lines of code to save the plot as a pdf
#dev.off()

### Results, Figures, and Tables from Appendix

## Reproducing Results from Table 3 in Appendix:
## Main Models of Trustworthiness of Mauritania/Burkina Faso

# Base Model

basemodel <- lm(trust_total ~ similaritycond + costcond + country
         + log(age) + male + education + vax + familyhadcovid
         + demdummy + repdummy + polorient
         + CIfa + MIfa + ISOfa + recipcond,
         data = d1)
summary(basemodel)

# Model 2

model2 <- lm(trust_total ~ similaritycond + costcond + country
          + log(age) + male + education + vax + familyhadcovid
          + demdummy + repdummy + polorient
          + ci + iso + mi + recipcond,
          data = d1)
summary(model2) 

# Model 3

model3 <- lm(trust_total ~ similaritycond + costcond + country
          + log(age) + male + education + vax + familyhadcovid
          + partyscale + polorient
          + mi + ci + iso + recipcond,
          data = d1)
summary(model3)

# For printing tables: 

stargazer(basemodel, model2, model3, 
          column.labels = c("Model 1", "Model 2", "Model 3"), 
          intercept.bottom = FALSE, 
          single.row=TRUE,  
          notes.append = FALSE,   
          star.cutoffs = c(0.05, 0.01, 0.001),  
          header=FALSE)


## Reproducing Results from Table 4 in Appendix:
## Interaction Models

# Reciprocity x Similarity 

recxsimmodel <- lm(trust_total ~ costcond + country
         + log(age) + male + education + vax + familyhadcovid
         + demdummy + repdummy + polorient
         + CIfa + MIfa + ISOfa + recipcond*similaritycond,
         data = d1)
summary(recxsimmodel)

# Reciprocity x Cost

recxcostmodel <- lm(trust_total ~ similaritycond + country
             + log(age) + male + education + vax + familyhadcovid
             + demdummy + repdummy + polorient
             + CIfa + MIfa + ISOfa + recipcond*costcond,
             data = d1)
summary(recxcostmodel)

# For printing tables: 

stargazer(recxsimmodel, recxcostmodel, 
          intercept.bottom = FALSE, 
          no.space=TRUE, 
          single.row=TRUE,
          notes.append = FALSE,   
          star.cutoffs = c(0.05, 0.01, 0.001),  
          header=FALSE,
          digits = 2) 

## Reproducing Results from Table 5 in Appendix:
## Specific Positive and Negative Reciprocity Expectations

# Positive Reciprocity Expectations

positivemodel <- lm(posrecip_total ~ similaritycond + costcond + country
         + log(age) + male + education + vax + familyhadcovid
         + demdummy + repdummy + polorient
         + CIfa + MIfa + ISOfa + recipcond,
         data = d1)
summary(positivemodel)

# Negative Reciprocity Expectations

negativemodel <- lm(negnecip_total ~ similaritycond + costcond + country
         + log(age) + male + education + vax + familyhadcovid
         + demdummy + repdummy + polorient
         + CIfa + MIfa + ISOfa + recipcond,
         data = d1)
summary(negativemodel)

# For printing tables: 

stargazer(positivemodel, negativemodel, 
          intercept.bottom = FALSE, 
          no.space=TRUE, 
          single.row=TRUE,
          notes.append = FALSE,   
          star.cutoffs = c(0.05, 0.01, 0.001),  
          header=FALSE,
          digits = 2) 

## Reproducing Figure 4: Familiarity with Countries in Our Sample

# Put data into new dataframe with long format, to make combining plots easier
d_long <- d1 %>%
  pivot_longer(
    cols = c(bffamiliar, mauritfamiliar), 
    names_to = "Category", 
    values_to = "Familiarity"
  )

#Create density plot using ggplot
familiarityplot <- ggplot(d_long, aes(x = Familiarity)) +
  geom_histogram(aes(y = ..density..),
                 bins = 11, colour = "black", fill = "white") +
  geom_density(aes(fill = Category), alpha = 0.3) + 
  scale_x_continuous(breaks = seq(0, 10)) +
  scale_fill_manual(
    values = c("bffamiliar" = "#35B779FF", "mauritfamiliar" = "#ED6925FF"),
    labels = c("bffamiliar" = "Burkina Faso", 
               "mauritfamiliar" = "Mauritania")
  ) +
  labs(x = "Familiarity", y = "Density", fill = "Country") +
  theme_classic() +
  facet_wrap(~Category, ncol = 2, 
             labeller = as_labeller(c(
               bffamiliar = "Burkina Faso",
               mauritfamiliar = "Mauritania"
             )))
#ggsave("familiarityplot.pdf", plot = familiarityplot, width = 16, height = 10, units = "cm")
#dev.off()

## Reproducing Results from Table 6 in Appendix:
## Effect of Familiarity on Trust

# Reciprocity x Familiarity with Burkina Faso

recxbf <- lm(trust_total ~ similaritycond + costcond + country
          + log(age) + male + education + vax + familyhadcovid
          + demdummy + repdummy + polorient
          + CIfa + MIfa + ISOfa + recipcond*bffamiliar,
          data = d1)
summary(recxbf)

# Reciprocity x Familiarity with Mauritania

recxmaur <- lm(trust_total ~ similaritycond + costcond + country
          + log(age) + male + education + vax + familyhadcovid
          + demdummy + repdummy + polorient
          + CIfa + MIfa + ISOfa + recipcond*mauritfamiliar,
          data = d1)
summary(recxmaur)

stargazer(recxbf, recxmaur, 
          intercept.bottom = FALSE, 
          no.space=TRUE, 
          single.row=TRUE,
          notes.append = FALSE,   
          star.cutoffs = c(0.05, 0.01, 0.001),  
          header=FALSE,
          digits = 2) 

